home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / randist / gamma.c < prev    next >
Encoding:
C/C++ Source or Header  |  2001-05-14  |  3.9 KB  |  167 lines
  1. /* randist/gamma.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <math.h>
  22. #include <gsl/gsl_math.h>
  23. #include <gsl/gsl_sf_gamma.h>
  24. #include <gsl/gsl_rng.h>
  25. #include <gsl/gsl_randist.h>
  26.  
  27. static double gamma_large (const gsl_rng * r, const double a);
  28. static double gamma_frac (const gsl_rng * r, const double a);
  29.  
  30. /* The Gamma distribution of order a>0 is defined by:
  31.  
  32.    p(x) dx = {1 / \Gamma(a) b^a } x^{a-1} e^{-x/b} dx
  33.  
  34.    for x>0.  If X and Y are independent gamma-distributed random
  35.    variables of order a1 and a2 with the same scale parameter b, then
  36.    X+Y has gamma distribution of order a1+a2.
  37.  
  38.    The algorithms below are from Knuth, vol 2, 2nd ed, p. 129. */
  39.  
  40. double
  41. gsl_ran_gamma (const gsl_rng * r, const double a, const double b)
  42. {
  43.   /* assume a > 0 */
  44.   unsigned int na = floor (a);
  45.  
  46.   if (a == na)
  47.     {
  48.       return b * gsl_ran_gamma_int (r, na);
  49.     }
  50.   else if (na == 0)
  51.     {
  52.       return b * gamma_frac (r, a);
  53.     }
  54.   else
  55.     {
  56.       return b * (gsl_ran_gamma_int (r, na) + gamma_frac (r, a - na)) ;
  57.     }
  58. }
  59.  
  60. double
  61. gsl_ran_gamma_int (const gsl_rng * r, const unsigned int a)
  62. {
  63.   if (a < 12)
  64.     {
  65.       unsigned int i;
  66.       double prod = 1;
  67.  
  68.       for (i = 0; i < a; i++)
  69.     {
  70.       prod *= gsl_rng_uniform_pos (r);
  71.     }
  72.  
  73.       /* Note: for 12 iterations we are safe against underflow, since
  74.      the smallest positive random number is O(2^-32). This means
  75.      the smallest possible product is 2^(-12*32) = 10^-116 which
  76.      is within the range of double precision. */
  77.  
  78.       return -log (prod);
  79.     }
  80.   else
  81.     {
  82.       return gamma_large (r, (double) a);
  83.     }
  84. }
  85.  
  86. static double
  87. gamma_large (const gsl_rng * r, const double a)
  88. {
  89.   /* Works only if a > 1, and is most efficient if a is large
  90.  
  91.      This algorithm, reported in Knuth, is attributed to Ahrens.  A
  92.      faster one, we are told, can be found in: J. H. Ahrens and
  93.      U. Dieter, Computing 12 (1974) 223-246.  */
  94.  
  95.   double sqa, x, y, v;
  96.   sqa = sqrt (2 * a - 1);
  97.   do
  98.     {
  99.       do
  100.     {
  101.       y = tan (M_PI * gsl_rng_uniform (r));
  102.       x = sqa * y + a - 1;
  103.     }
  104.       while (x <= 0);
  105.       v = gsl_rng_uniform (r);
  106.     }
  107.   while (v > (1 + y * y) * exp ((a - 1) * log (x / (a - 1)) - sqa * y));
  108.  
  109.   return x;
  110. }
  111.  
  112. static double
  113. gamma_frac (const gsl_rng * r, const double a)
  114. {
  115.   /* This is exercise 16 from Knuth; see page 135, and the solution is
  116.      on page 551.  */
  117.  
  118.   double p, q, x, u, v;
  119.   p = M_E / (a + M_E);
  120.   do
  121.     {
  122.       u = gsl_rng_uniform (r);
  123.       v = gsl_rng_uniform_pos (r);
  124.  
  125.       if (u < p)
  126.     {
  127.       x = exp ((1 / a) * log (v));
  128.       q = exp (-x);
  129.     }
  130.       else
  131.     {
  132.       x = 1 - log (v);
  133.       q = exp ((a - 1) * log (x));
  134.     }
  135.     }
  136.   while (gsl_rng_uniform (r) >= q);
  137.  
  138.   return x;
  139. }
  140.  
  141. double
  142. gsl_ran_gamma_pdf (const double x, const double a, const double b)
  143. {
  144.   if (x < 0)
  145.     {
  146.       return 0 ;
  147.     }
  148.   else if (x == 0)
  149.     {
  150.       if (a == 1)
  151.     return 1/b ;
  152.       else
  153.     return 0 ;
  154.     }
  155.   else if (a == 1)
  156.     {
  157.       return exp(-x/b)/b ;
  158.     }
  159.   else 
  160.     {
  161.       double p;
  162.       double lngamma = gsl_sf_lngamma (a);
  163.       p = exp ((a - 1) * log (x/b) - x/b - lngamma)/b;
  164.       return p;
  165.     }
  166. }
  167.